Rare Transition Events in Nonequilibrium Systems with State-Dependent Noise: 
Application to Stochastic Current Switching in Semiconductor Superlattices 



O 

(N 
-i— > 
O 

O 



> 

m 
O 

od 
O 
o 



Matthias Heymann, 1 '* Stephen W. Teitsworth, 2, t and Jonathan C. Mattingly 1 '* 

1 Duke University Mathematics Department, Durham, North Carolina 27708, USA 
2 Duke University Physics Department, Durham, North Carolina 27708, USA 
(Dated: October 19, 2010) 

Using recent mathematical advances, a geometric approach to rare noise-driven transition events in 
nonequilibrium systems is given, and an algorithm for computing the maximum likelihood transition 
curve is generalized to the case of state-dependent noise. It is applied to a model of electronic 
transport in semiconductor superlattices to investigate transitions between metastable electric field 
distributions. When the applied voltage V is varied near a saddle-node bifurcation at Vth, the mean 
life time (T) of the initial metastable state is shown to scale like log(T) oc \Vth — V\ 3 ^ 2 as F / Vth- 
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The study of noise-driven transitions between 
metastable states is of current interest for a range of sys- 
tems, for example, chemical reaction systems [1], nano- 
and micromechanical oscillators [2-4], magnetic tunnel 
junctions [5], and biochemical networks [6]. Accordingly, 
the development of mathematical and numerical methods 
[7-10] for understanding in particular the mean transi- 
tion (or escape) time and the most probable escape path 
is an active area with implications for a variety of fields. 
Since the widely-used string method [7] is restricted to 
the case of gradient drift, the geometric minimum action 
method (gMAM) [9-11] was recently introduced for non- 
gradient systems (as a successor of the minimum action 
method, MAM [8]). However, none of its implementa- 
tions can at present deal with multiplicative (i.e., state- 
dependent) noise, leaving it inapplicable to a number of 
important physical problems. 

In this paper we fill this gap by extending gMAM 
to systems with multiplicative noise, and we also intro- 
duce a useful new method for locating unstable equilib- 
rium points in high-dimensional systems. The mathe- 
matical framework, large deviation theory [12], is pre- 
sented in a non-standard geometric way (i.e., based on 
zmparameterized curves), which allows us to incorporate 
a key result of a recent study [13] proving the existence 
of a maximum likelihood transition curve 7*. 

These techniques are illustrated by applying them to 
the problem of stochastic current switching in a semicon- 
ductor superlattice (SL), which is a linear array of N > 1 
quantum wells that arc coupled to each other via elec- 
tron tunneling. While the deterministic nonlinear prop- 
erties of SLs have been studied extensively and are gener- 
ally well- understood [14], very little theoretical work [15] 
exists for explaining experimental measurements [16] of 
noise-induced current switching in this system. One ob- 
ject of interest is the mean escape time (T) from the 
metastable states in the SL system, and its dependence 
on key control parameters such as the applied voltage 
V. Since the noise in the SL model is multiplicative, 
previous studies could not yet use gMAM and instead 



relied on direct simulation of the underlying stochastic 
differential equations (SDEs) [15]. The resulting lack of 
accuracy and high computational cost (note that (T) be- 
comes exponentially large in the zero-noise limit) can now 
be circumvented by the non-trivial extension of gMAM 
introduced here. 

This allows us, for the first time, to determine the scal- 
ing behavior of (T) for the SL system. Specifically, we 
show that log(T) cx \V th - V\ 3 ^ 2 as V / V th , where V th is 
the value above which the initial metastable state ceases 
to exist. Previously, this scaling result was only estab- 
lished for the one-dimensional case, i.e., for the case of a 
single quantum well with small cross-sectional size [17]. 

We begin by presenting the SDE model for the SL sys- 
tem. Electronic transport properties of weakly-coupled 
superlattices are accurately described by a spatially dis- 
crete model that incorporates sequential resonant tun- 
neling between successive quantum wells [14, 18]. For a 
superlattice with N quantum wells, this model consists 
of the Poisson and the charge continuity equations 
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where Fi denotes the spatially-averaged electric field in 
the i th period of the SL, m is the 2D electron density in 
the i th well, e is the elementary charge, e is the dielectric 
constant, Nd = 1.5 x 10 11 cm -2 is the 2D doping density, 
and Ji is the tunneling current density from the i th to 
the (i + l) th well [19]. Differentiating (1) with respect to 
time and inserting the result into (2) gives 
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where the total current density J(t) is the same for all 
periods. The bias condition is 
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FIG. 1. (a) The /(Fi/F max ) curve, (b) The J-V curve. Every 
point on the two solid branches corresponds to an attractor 
of our system, points on the dotted branch correspond to a 
saddle point. The vertical line at Vth marks the location of 
the saddle node bifurcation. 



where V is the time-independent voltage bias across the 
entire SL (which is the control parameter for this system), 
and where I is the distance between successive wells. 

The current density Jj is developed using a transfer 
Hamiltonian, and it can be written as a function of the 
local field and the neighboring charge densities [18]: 
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for i = 1, . . . , N — 1, where c\ = 1.68 • 10 10 cm 2 , C2 = 
3.0145 cm/kV, v M = 1.691 m/s, F max = 3.945 kV/cm, 
and where the function f(Fi/F mElx ) is shown in Fig. 1(a) 
[20]. For the present study, we use standard Ohmic 
boundary conditions 



Jo = gF and J N = gF N 
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where g = 0.08 (51m) -1 denotes the effective contact 
conductivity. 

Equations (1) and (3)- (6) constitute a complete dy- 
namical system in (Fq, . . . , Fn) describing electronic 
transport in the SL system and, for the parameters used 
here, are well-known to yield two stable points within 
a range of bias values V [18]. By (3), (5)-(6) and (1), 
for each equilibrium solution we have Jo = • • • = Jn = 
J = const, and so the bistability can be illustrated 
by the current- voltage curve in Fig. 1(b) which for ev- 
ery equilibrium point shows the associated value for J. 



The full curve consists of N + 1 stable current branches, 
the first N of which each ends in a saddle node bifur- 
cation; Fig. 1(b) shows only the 4 th and 5 th branches 
(the two solid curves), which are the focus of our study. 
Within the n th current branch (counted from V = 0), 
every stable state (Fq, . . . , Fn) consists of two domains 
(Fq, . . . , Fn-u) and (Fn-u+i, ■ ■ ■ , Fn) within which the 
entries Fi are roughly constant (low in the first and high 
in the second domain). 

Let us now consider the effect of shot noise, which re- 
sults from the discrete and random nature of the tun- 
neling electrons [15, 21]. This noise induces rare tran- 
sitions from the metastable state on branch n to the 
one on branch n + 1, which causes the overall current 
to jump and the high field domain to expand by one pe- 
riod. Adding this noise to the tunneling current density, 
we modify (3) as 



J(t)=e^ + J i + Jf\t), 



0,...,N. (7) 



Using that the noise terms J\ (t) are delta-correlated in 
time and space, one can solve (7) for dFi to write [21] 
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for i = 0,...,N, where Bf = ^ (here a denotes the 
cross-sectional area of the SL), and where the W % are 
independent Brownian motions. 

In order to apply gMAM to this model, we need to 
resolve the problem that the noise terms in the (N + 1)- 
dimensional system (8) are dependent (note that, be- 
cause of (4), summing the second bracket in (8) over 
all i gives 0). We will therefore instead consider the N- 
dimensional transformed system for the variables n,;, 

dn, = - (dFi - dFi_x) 
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for i = 1, . . . , N. To see that this is a closed system, 
observe that the linear system for the Fi given by (1) 
and (4) is solved by 
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for i = 0, . . . , N, and so by (5)-(6) we actually have 
Ji = Ji{ni, ■ ■ ■ ,it>n) for Vi = 0, . . . , N. Denoting X := 
(m, . . . , tin), the system (9) is thus of the form 

dX? = b(X?) dt + V^er(A7) dW t , (10) 
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where iWt)t>o is an V-dimensional Brownian motion, 
77= (ea) _1 is a small parameter, and where the drift vec- 
tor field b(X) and the diffusion matrix a{X) are given by 



(ii) 



(12) 



b(X) = -(Jo — Ji, . . . , Jiv-i — Jjv) 
e 

A/Jo ••• 

V ■•• v/j^TT -v^y 



with Ji = Ji(X). As we will see below, this solves our 
problem of degenerate noise. 

One of the central results of large deviation theory is 
that the mean time for a transition from one metastable 
state x\ to another, X2, in a system of the form (10), is 
given by (T v ) « e u{x u x 2 )/ji as 77 X^ 0, where U is the 
quasipotential [12]. It is well-known that in the gradient 
case (i.e., if b — —VU for some potential U, and if a is 
the identity matrix) we have U{xx,X2) = U(x s ) — U(xi), 
where x s is the lowest-energy saddle point between the 
two basins of attraction; the maximum likelihood tran- 
sition path is in this case the path of minimum en- 
ergy, which can be numerically computed using the string 
method [7]. However, the theory extends to the general 
case (10), and the quasipotential U is then given by 

U(xi,x 2 ) = inf / (\b(x)\ x \dx\ x - (b(x),dx) x ), (13) 

where T^ 2 is the set of rectifiable (i.e., finite- length) 
curves leading from x\ to X2, and where (u,v) x := 
(u, A(x)~ 1 v) and \u\ x := (u,u)^ 2 for Vw, v e R N . Here 
we define A{x) := a(x)a(x) T , which we require to be 
positive definite for every x. 

Equation (13) is a geometric reformulation [9] of a 
more common formula for U that is based on time- 
parameterized paths [12, Eqns. (3.5) and (4.1)], and it 
has gotten some attention recently for its analytical and 
numerical advantages. In particular, the infimum in (13) 
is typically achieved at some minimizing curve 7* G Y^l 
[13], i.e., the maximum likelihood transition curve, and 
so denoting the integral in (13) by £(7), we have 



U{x u x 2 ) =S( 7 *). 



(14) 



Furthermore, (13) is the basis of an efficient algorithm 
(the geometric Minimum Action Method, or gMAM) for 
computing 7*. The algorithm works by moving some 
initial curve 70 € T®* successively into the direction of 
steepest descent while keeping the end points fixed. If 
we write 

5( 7 ) = S(<p) = [ (|%)|>% - {b(<p),<p') v ) da 



for some parameterization ip(a) : [0, 1] — > R N of 7 then 
gMAM numerically solves the PDE d T ip(T, a) = —SS(ip), 
(t, a) £ [0, 00) x [0, 1], subject to the boundary conditions 
ip(r = 0, • ) = ipo, (p(-,a = 0) = xt, tp(-,a = l) = x 2 . 
After preconditioning [9], this PDE can be written as 

d T tp = \ 2 <p" - X(Wb{ip) + C)cp' 

+ A{tp)(Vb{<p) + ±C) T 8 + \\'^, (15) 

where ip' and ip" denote the a-derivatives of tp, X((p,ip') 
:= \b(cp)U\<p%, Ofay?) := A{<p)-\\<f/ - %)), and 
C((p,ip') is the (N x iV)-matrix whose i th column is 
[c^A] 1^ _ v 0. To stabilize the algorithm, a semi-implicit 
code is used for the if" term, and the discretization points 
are redistributed equidistantly along the curve after each 
iteration step. See [9, 10] for the derivation of (15) in 
the additive case (i.e., if A is the identity matrix) and for 
further details on the implementation of gMAM. 

To see that gMAM is applicable to our system (10)- 
(12), first note that 

fJa + Jx -Ji ■•• \ 
— Ji J\ + J 2 —J2 '■. : 
A = aa T = . . ' ■ . 

: '• — Jjv-i 

\ ••• -J N -i Jn-i + JnJ 

it is then straightforward to show by induction on N 
that det(^4) = X^fclo YliLo i=tk which is positive since 
Ji > for Vi = 0, . . . , N. 

To apply gMAM to the SL system, we first fix V = 
0.52 Volt (see Fig. 1(b)) and locate the two stable points 
(m,...,niv) of the transformed system with Newton's 
method. To find the saddle point in between, we make 
a novel application of the string method [7]. Briefly, the 
string method is a fast and easy algorithm for finding a 
curve that connects the two attractors and is everywhere 
parallel or antiparallel to the drift b. While the original 
purpose of this algorithm was to find transition curves in 
gradient systems (which have this geometric property), 
in general systems one can still make use of the fact that 
by construction such a curve must contain an equilib- 
rium point. One can therefore find the saddle point by 
(i) applying the string method, (ii) searching for an equi- 
librium point along the obtained curve, and (iii) increase 
its accuracy by using Newton's method. 

Finally, we find the minimum action curve 7* leading 
from the first attractor to the saddle point, and we com- 
pute the associated action 5(7*). (Note that, as a con- 
sequence of the non-gradient nature of this system, the 
curve 7* is found to deviate significantly from the unsta- 
ble manifold of the saddle point, see also [4].) By increas- 
ing V and using continuity, we determine the three equi- 
librium points, and then 7* and S ("/*), all as functions 
of V. We identify V th = 0.5578288 Volt as the value V 
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FIG. 2. (a) S(V), (b) S' C (V). The solid circles are the values 
found by gMAM, solid and dotted curves are the approxima- 
tions (16) and (17), respectively. The vertical lines mark Vth- 



at which the attractor on the upper branch collides with 
the saddle point . 

Fig. 2(a) shows plots of the action S(V) = S(-y*(V)). 
The values found by gMAM are shown as solid circles, 
the solid and dotted curves are the approximations found 
below. As the log-log plot in the inset illustrates, near the 
threshold Vth we obtain the lowest-order approximation 



S(V)^s v p , where v = \V th - V\/V th , 



(16) 



/3 = 1.5±0.0001 and s = (1.905±0.0002)-10 30 C^cm" 2 . 
The arrows mark the data point at which the slope of the 
log-log plot increases to 1.6, and so we find that in this 
sense (16) is valid for 0.55 Volt < V< V th . 

To obtain higher-order terms, we consider S C (V) := 
S(V)/(s u 3 / 2 ), which by construction fulfills S c (V th ) = 1. 
The linear approximation of S' C (V) in Fig. 2(b) leads to 
a quadratic approximation of S C (V), which yields 



S(V) 



SoV 3/2 + SiW 5/2 



S 2 V 



7/2 



(17) 



where Sl = 2.35 • 10 3i C^cm" 2 and s 2 = 10 32 C^cm" 2 . 
This appears to be a good fit in the entire voltage range 
0.52 Volt < V< V th , see the dotted line in Fig. 2(a). 

In this paper, we have generalized the gMAM algo- 
rithm for finding transition curves and their associated 
actions to SDE systems with multiplicative noise. We 
have used a new method for locating saddle points, and 



we have demonstrated the use of a geometric formulation 
(13) of the quasipotential. 

We have applied our techniques to the problem of 
stochastic current switching in semiconductor SLs, and 
to do so, we transformed the standard SL model into 
a form with non-degenerate noise matrix. Our results 
for the SL show that in the zero-noise limit we have 
log(T) oc \V th - V| 3/2 as V /■ V th , and we are able to 
quantify the range in which this scaling is valid. In ad- 
dition, we have found a higher-order approximation that 
is valid also for larger values of \Vth — V\. 
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